00100	C This is taken from the "Handbook of Mathematical Functions"
00200	C by Abramowitz and Stegun, Dover press S1272, 1965.
00300	C Chapter 9 is on Bessel Functions. The polynomial
00400	C approximations are found in section 9.4, pages 369
00500	C and 370. The recursion formula for generating higher
00600	C orders from J0 and J1 is found at the bottom of
00700	C table 9.1, page 390
00800	
00900		SUBROUTINE BESCOEF(MI,X,I,K)
01000		REAL MI,J0,J1
01100		COMMON FREQ(3,0/50,50),FUNC(50)
01200		IF(I.GT.1)GO TO 30
01300		IF(MI.GE.3.0)GO TO 10
01400		xs=(MI/3)**2
01500		j0=1.0+xs*(-2.2499997+xs*(1.2656208+xs*(-0.3163866+
01600		1xs*(0.0444479+xs*(-0.0039444+xs*0.00021)))))
01700		j1=MI*(.5+xs*(-0.56249985+xs*(0.21093573+
01800		1xs*(-0.03954289+xs*(0.00443319+xs*(-0.00031761+
01900		2xs*0.00001109))))))
02000		GO TO 20
02100	10	xs=3.0/MI
02200		j0=(1.0/sqrt(MI))*(0.79788456+xs*(-0.00000077+
02300		1xs*(-0.0055274+xs*(-0.00009512+xs*(0.00137237+
02400		2xs*(-0.00072805+xs*0.00014476))))))*
02500		3cos(MI-0.78539816+xs*(-0.04166397+
02600		4xs*(-0.00003954+xs*(0.00262573+
02700		5xs*(-0.00054125+xs*(-0.00029333+
02800		6xs*0.00013558))))))
02900		j1=(1.0/sqrt(MI))*(0.79788456+xs*(0.00000156+
03000		1xs*(0.01659667+
03100		1xs*(0.00017105+xs*(-0.00249511+xs*(0.00113653-
03200		2xs*0.00020033))))))*
03300	 	3cos(MI-2.35619449+xs*(0.12499612+xs*(0.0000565+
03400		4xs*(-0.00637879+xs*(0.00074348+xs*(0.00079824-
03500		5xs*0.00029166))))))
03600	20	FREQ(2,0,K)=J0
03700		FREQ(2,1,K)=J1
03800		FREQ(2,2,K)=J1
03900		RETURN
04000	30	Y=I
04100		IF(MI.LE.0.0000001)GO TO 50
04200		IJ=I*2		
04300		X=((2.0*(Y-1.0))/MI)*FREQ(2,IJ-2,K)-FREQ(2,IJ-4,K)
04400		RETURN
04500	50	X=0
04600		RETURN
04700		END
     

00050		SUBROUTINE SIDEBANDS(CAR,XMOD,CINDEX,K)
00100		COMMON FREQ(3,0/50,50),FUNC(50)
00150		IT=-1.0
00300		MAX=CINDEX+5.0
00350		IF(FREQ(1,50,1).LT.MAX*2.)FREQ(1,50,1)=MAX*2.
00400		DO 100 J=0,MAX*2,2
00500		I=J/2
00600		NORDER=I
00700		IF(I.NE.1)CALL BESCOEF(CINDEX,COEF,NORDER,K)
00800		IF(I.GT.0)GO TO 20
00900		FREQ(1,I,K)=CAR
01000		FREQ(3,I,K)=0.0
01100		GO TO 30
01200	20	XI=I
01300		YMOD=XI*XMOD
01400		FREQ(1,J-1,K)=CAR+YMOD
01500		FREQ(1,J,K)=CAR-YMOD
01600		IF(I.EQ.1)GO TO 10
01700		FREQ(2,J-1,K)=COEF
01800		FREQ(2,J,K)=COEF
01900	10	FREQ(3,J-1,K)=I
01950		IF(IT)I=-I
02000		FREQ(3,J,K)=I
02050		IT=-IT
02100	30	CONTINUE
02200	100	CONTINUE
02500	102	FORMAT(3F)
02525		NAX=MAX
02537		CALL REFLECT(NAX,K)
02550		RETURN
02600		END
02700	
     

00100		SUBROUTINE REFLECT(MAX,K)
00200		COMMON FREQ(3,0/50,50),FUNC(50)
00210		DO 100 N=0,MAX*2
00255	100	IF(FREQ(3,N,K).LT.0.0)FREQ(2,N,K)=-FREQ(2,N,K)
00300		DO 201 J=0,(MAX*2)-1
00400		DO 201 I=J+1,MAX*2
00500		IF(FREQ(1,J,K).OR.FREQ(1,I,K).EQ.99999.)GO TO 201
00600		IF(ABS(FREQ(1,J,K)).NE.ABS(FREQ(1,I,K)))GO TO 201
00700		IF(FREQ(1,J,K).GT.FREQ(1,I,K))GO TO 20
00800		FREQ(2,I,K)=FREQ(2,I,K)-FREQ(2,J,K)
00900		FREQ(1,J,K)=99999.
01000		GO TO 201
01100	20	FREQ(2,J,K)=FREQ(2,J,K)-FREQ(2,I,K)
01200		FREQ(1,I,K)=99999.
01300	201	CONTINUE
01400		RETURN
01500		END
01600	
01700	
01800	C     	MAIN PROGRAM
01900		COMMON FREQ(3,0/50,50),FUNC(50)
02100	10	ACCEPT 100,C,XM,ZI1,ZI2
02200	100	FORMAT(4F)
02250		CALL GEN
02275		DO 300 IS=1,25
02287		CI=ZI1+(ZI2-ZI1)*FUNC(IS)
02300	300	CALL SIDEBANDS(C,XM,CI,IS)
02400		CALL DISP
02500		MAX=FREQ(1,50,1)
03410	C	TYPE 101,(((FREQ(K,L,M),K=1,3),L=0,MAX),M=1,25)
03430	101	FORMAT(1H (3F/))
03450		GO TO 10
03500		END
     

03600	
03700		SUBROUTINE GEN
03800		COMMON FREQ(3,0/50,50),FUNC(50)
03900		X=1./25.
04000		DO 100 I=1,25
04100		Y=I-1
04200	100	FUNC(I)=X*Y
04300		RETURN
04400		END
04500	
04550		SUBROUTINE DISP
04555		DIMENSION I(1),IJ(1000)
04560		COMMON FREQ(3,0/50,50),FUNC(50)
04600		CALL DPYSET(1,IJ,1000)
04700		CALL ALINE(-400,0,400,0)
04800		CALL ALINE(-400,400,-400,0)
04900		MAX=FREQ(1,50,1)
04950		KL=1
05000		DO 200 J=0,MAX
05005		IF(FREQ(1,J,KL).EQ.99999.)GO TO 200
05010		IX=ABS(FREQ(1,J,KL))-400.
05012		ZZ=IX
05015		IY=(ZZ+400.)*(-1.)+400.*FREQ(2,J,KL)
05017		BASE=(ZZ+400.)*(-1.)
05020		KL=KL+1
05022		CALL AIVECT(IX,IY)
05025		DO 200 K=2,25
05100		IF(FREQ(1,J,K).EQ.99999.)GO TO 200
05102		IF(FREQ(1,J,K).EQ.0.0)GO TO 200
05105		IX=IX+20
05110		IY=FREQ(2,J,K)*400.+BASE
05112		IBASE=BASE
05114		IF(IY.LE.IBASE)IY=(IABS(IY)-IABS(IBASE))+IBASE
05115		CALL AVECT(IX,IY)
05600	200	CONTINUE
05700		CALL DPYOUT(1)
06000		RETURN
06100		END